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ABSTRACT 

Pulsar timing observations have revealed companions to neutron stars that include other neutron 
stars, white dwarfs, main-sequence stars, and planets. We demonstrate that the correlated and ap- 
parently stochastic residual times of arrival from the millisecond pulsar B1937+21 are consistent with 
the signature of an asteroid belt having a total mass < 0.05 M®. Unlike the solar system's asteroid 
belt, the best fit pulsar asteroid belt extends over a wide range of radii, consistent with the absence 
of any shepherding companions. We suggest that any pulsar that has undergone accretion-driven 
spin-up and subsequently evaporated its companion may harbor orbiting asteroid mass objects. The 
resulting timing variations may fundamentally limit the timing precision of some of the other mil- 
lisecond pulsars. Observational tests of the asteroid belt model include identifying periodicities from 
individual asteroids, which are difficult; testing for statistical stationarity that become possible when 
observations are conducted over a longer observing span; and searching for reflected radio emission. 
Subject headings: planets — pulsars: general — pulsars: specific (PSR B1937+21) — stars: neutron 



1. INTRODUCTION 

Arrival time measurements of rotation-powered pulsars 
show time-correlated excesses after accounting for con- 
tributions from the spin-down of the p ulsar, astrometric 
varia tions, and obvious orbital motion (jCordes fc Downs! 
Il985h . often referred to as timing noise (TN), but which 
we will refer to as red noise (RN). In slowly spinning 
canonical pulsars (CPs, pulsars with spin frequencies 
0-5 ^ v ;$ 30 Hz, and relatively high surface magnetic 
fields B ~ 10 12 G), the excess is likely caused by ro- 
tational irregularities associated wi th a combination of 
magnetospheric t orque fluctuations (|Kramer etaIll20M 
iLvne et"al1l2010h and instabilities arising from differen- 
tial rotation of the neutron star (NS) crust and core 
(|Jonedll99Ch . 

Pulsars that have undergone accretion-driven spin- 
up and magnetic field attenuation are much more sta- 
ble. The stability of these millisecond pulsars (MSPs, 
30 < v i5 700 Hz) has led to the discovery of the 
first Earth-mass planets outside of the solar system 
(Wols zczan fc Fraill IT992T ) . stringent tests of general rel- 
ativity through the study of g ravitational intera ctions 
with white dwarf companions (jFreire et all I2012f ). and 
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const raints on nuclear equations of state (jDemorest et al.1 
l201Ch . If timing precision can be improved further it may 
be possible to detect the signature of gravitational waves 
passing through the solar neighborhood using an ensem- 
ble o f pulsars (i.e., pulsar timing array. iFoster fc Backer! 
fl990lh 

In some pulsars, the timing residuals contain tempo- 
rally correlated stochastic noise that often appears con- 
sistent with a 'red' power spectrum that has significant 
power at low fluctuation frequencies. The most striking 
example of this in an MSP is seen in the timing analysis of 
PSR B1937+21. In Figure [Q the post-fit residual times 
of arrival (TOAs) spanning 26 yr are displayed for this 
object. The observations used to form these residuals 
came from the Arecibo, Effelsberg, Nangay, and Wester- 
bork telescopes. The observations and the model used 
to derive the residual TOAs are further described in Ap- 
pendix [A] The residuals show correlated variations that 
are similar in spectral nature to that observed in CPs, 
but at a much smaller level. The root mean square (rms) 
strength of the RN for the observations presented in Fig- 
ure Q] is <trn,2 ~ 45 //s, where the subscript '2' implies 
the rms residual after a second-order polynomial fit. In 
contrast, the rms strength of RN observed in many CPs 
can exceed 1 s over a similar observing span. 

This excess can be interpreted in a number of ways. 
The residuals could be associated with rotational insta- 
bilities similar to those that cause RN in CPs that we 
will refer to as red noise. Indeed, a self consistent scaling 
relationship relating RN to spin frequency and frequency 
derivative has been developed that explains the measured 
levels and up per limits of RN in MSP s and the slower 
spinning CPs JSh annon fc Cordesll2010l ). which suggests 
that RN may be associated with the spin properties of 
the NS. 

Additional contributions arise from refraction and 
scattering of radio waves in the interstellar medium 
(ISM). However, these effects are highly wavelength- 
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dependent. After correction the timing residuals of 
PSR B1937+21 show a persistent achromatic compo- 
nent (|Cordes et al.l [1990T) . that we verify here. Achro- 
matic effects arise from uncertainties in the solar sys- 
tem ephemeris used to tran sform topocentric TQA s to 
the solar system barycenter ([Champion et al.ll2010( ) and 
from t he perturbation s from long-period gravitational 
waves ([Detweilerlll979D . but these effects are small and 
are correlated between pulsars. The residuals seen in 
PSR B1937+21 are much larger than those seen in other 
MSPs so these effects are secondary if not negligible. 
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Figure 1. Residual TOAs At for PSR B1937+21 over the 26 yr 
observing span used in this analysis. 
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Figure 2. Maximum entropy power spectrum for the observa- 
tions of PSR B1937+21 presented in Figure □ The methods for 
generating power spectra are discussed in Appendix|B] The plateau 
in power at / < 0.1 yr - 1 is associated with quadratic fitting intro- 
duced in TOA modeling, which high-pass filters the time series. 



Here we investigate an alternative mechanism for the 
correlated microsecond timing residuals in MSPs: physi- 
cal displacement of the pulsar due to recoil from circum- 
pulsar objects. While the residuals from PSR B1937+21 
are too small and too incoherent to be associated with 
a few large planets in compact orbits, they could be as- 
sociated with a larger ensemble of asteroid-mass bodies 
that collectively produce the residuals. 

There is direct evidence that planetary and pro- 
toplanetary environments exist around some pul- 
sars. Most obviously, there is the planetary system 
around the MSP B1257+12 (jWolszczan fc Fraill Tl992t 
iKonacki fc Wolszczanl l200l comprising three moon to 
Earth-mass planets in ss 25 d, 66 d, and 98 day orbits 
that induce an rms scatter of k, 1 ms on the residual 
TOAs. There is also evidence for a planet in a wide orbit 
around the globular cluster pulsar B1620— 26 and its bi- 
nary white dwarf companion, inferred from both secular 
changes in the spin d own of the pulsar and evolution of 
the WD-pulsar orbit ([Backer et al.lll993t iThorsett etHI 
1999). This planet was likely captured in a three-body 
interaction comm on in the dense cluster environment 
(Sigurdssonl ll993T ) representing a formation channel that 
will not be discussed further here. 

There is evidence for disks and planetary systems 
around other neutron stars as well. The young 
magnetar 4U 0142+61 shows excess mid-infrared flux 
that can b e attributed to t hermal emission from a 
dust disk ([Wang et al.l 120061 ). However, the excess 
has also been attributed t o ma gnetospheric emission 
(|Beloborodov fc Thompson! 12007ft . More recently, the 
impact of an asteroid-mass object on a neutron star has 
been proposed as the cause of an unusual gamma-ray 
burst ([Campana et al.ll2011l ). 

Circumstantial evidence comes from the likely pres- 
ence of asteroid belts around another class of compact 
objects: white dwarf stars (WDs). Accretion of aster- 
oids is thought to pollute atmosphere s and cause anoma- 
lously high metallicity in some WDs (|Koester fc Wilkenl 
120061 : iFarihi et aTll2011h . Tidal disruption of asteroids 
is hypothesized to cause d ust disks around other WDs 
([Debes fc Sigurdssonll2002l ). While the planetary nebula 
phase that precedes the white dwarf end state for low 
mass stars is not as traumatic as the supernova explosion 
that precedes the neutron star end state for higher mass 
stars, the presence of asteroids around WDs suggests that 
rocky bodies can exist around post main sequence stars 
in relatively extreme environments. 

We report our findings as follows: in Section [21 we 
discuss the perturbations to pulsar TOAs caused by as- 
teroid belts, in Section [3l we discuss the formation of 
asteroidal disks around pulsars, in Section [4[ we show 
that many configurations of asteroid belts can be stable 
for the entire many Gyr lifetime of an MSP system; in 
Section O we demonstrate through simulation that there 
are many plausible disk configurations that explain the 
observed residual TOAs of PSR B1937+21; in Section 
O we compare the PSR B1937+21 system to other sites 
of planet formation; in Section we motivate further 
observations that can be used to assess the plausibility 
of the asteroid belt model; and in Section [51 we discuss 
implications for precision pulsar timing. 
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2. TIMING PERTURBATIONS FROM CIRCUMPULSAR 
OBJECTS 

The reflex motion of a pulsar due to a companion 
causes pulse TOAs to vary. This effect has been used 
to detect comp anion objects ranging from massive main 
seq uence stars dJohnston et al.ll992f) to Earth mass plan- 
ets (jWolszczan fc FraiW1992lP " 

In this section, we first estimate the rms TOA pertur- 
bation otoa produced by an orbiting body for compar- 
ison with the published rms residual, which enables us 
to make order-of-magnitude assessments of the size and 
structure of the system surrounding PSR B1937+21. In 
later sections, we will compare the power spectrum of the 
observed residuals to those induced by simulated asteroid 
belts. 

For a single object with a mass to orbiting a neutron 
star of mass Mns in a circular orbit with a radius a, the 
rm s residual TO A perturbation after a second order fit 
is (|Cordeslll993lfl 



<?2 ~ F(P mh ,T) 



a sin i \ / m 



V2 J \Mi 



NS 



0.85 ms F(P OTh ,T) 
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Mm 



1 yr 



2/3 



sini, (1) 



where i is the inclination angle of the orbit with respect 
to the line of sight. For a given mass, objects at larger 
radii produce larger timing perturbations. In the second 
line we have expressed the perturbation as a function of 
orbital period. 

The transmission function P(P or b,T) accounts for the 
TOA fitting procedure used to model pulsar spin phase, 
frequency and frequency derivative, which removes any 
parabolic perturbation to pulsar TOAs. In Figure [3] 
we show P(P or b,P) for our observing cadence, derived 
through simulations described in Appendix [C] Ignoring 
the astrometric terms (which absorb power in two nar- 
row windows at orbital periods of 1/2 and 1 yr, which 
are unimportant to this analysis), we find an empirical 
relation 



p(p ) (/Wl9yrr 2 - 83 
° lh> V 1 + (W19 yr)" 5 



67 ' 
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which is displayed as the dashed line in Figure [3] We find 
that when P or b <C 19 yr, F ~ 1 and when P or b 3> 19 yr, 



FocP 



-2.83 



oc a 



-4.25 



The effect of the outer regions on 



the disk is suppressed by the necessary quadratic fitting. 
Nonetheless, the outer regions can have significant influ- 
ence on the residual TOAs if the objects have sufficiently 
large mass. 

For PSR B1937+21, the lack of an obvious periodicity 
implies that there is no single, sufficiently massive com- 
panion with a period shorter than approximately 26 yr, 
which is the observing span of our observations. The 
total mass in a large ensemble of asteroids at random or- 
bital phases can exceed that of a single planet by several 
orders of magnitude because the residuals are an incoher- 
ent sum of the individual sinusoidal contributions. For 

9 Our expr ession differs by a factor of l/y/2 from that given in 
ICordesI 1119931) only because we present the rms residual rather than 
the peak-to-peak value, and we include the transmission factor F. 
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Figure 3. Transmission function as a function of P or b, and semi- 
major axis a, derived from a set simulations described in Appendix 
ICl The solid lines show the ±l-cr standard deviations about the 
mean (not displayed). The dashed line shows the empirical rela- 
tionship that we have adopted in Equation (T2J. 



j = l,N a asteroids, the rms perturbation to the TOAs 
is 



otoa 



N a 



1/2 



(3) 



If all asteroids have the same inclination, as would be 
the case if they originate from a disk, and we make the 
simplifying assumption that the masses and semi-major 
axes of the asteroids are statistically independent, the 
rms of the time series can be expressed as 



otoa 



where 
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100 
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(6) 



/a is the probability density function (PDF) for a, and 
Jm is the PDF for to,. 

By defining the total mass of the asteriod belt in terms 
of the mean asteroid mass, Mbeit = N a (m), we express 
the rms residual as 

(N a \- 1/2 f(a 2 F 2 ) 1 / 2 \ /M hclt 
otoa - 0.76 ms (sin, f — 1 ( ^AU - J [Tm^ 
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where £ is the ratio of the rms to mean mass: 



c 



(/?? 



2\l/2 
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If the mass probability density function follows a power- 
law, i.e. /m(^) oc m~ a , in the mass range M m \ n to 



> M n 



(2 -a) 



V(3-o)(l-a) 

(2 -a) 
V(3-a)(a-l) 



a < 1 



(a-l)/2 



1 < a < 2, 



C 9 ) 

The values of £ show marked variation depending 
on the configuration. If the disk has a relatively flat 
mass distribution or a narrow mass distribution with 
-Wmin ~ M max , then (« 1. A narrow mass distribution 
could arise if a dynamically cold disk of planetesimals 
was formed after a runaway growth phase, a scenario we 
discuss below. 

It is possible that the objects have a wide range of 
masses with a steep distribution. If the syst em is sup- 
porte d by collisional interactions, a — 11/6 (Dohnanvil 
1969) and objects could have a wide range of masses. 
Observations indicate that the size distribution of So- 
lar system objects has a distribution that is somewhat 
less steep with a ps 1.5 (jTerai fe Itobl l2011h . In both 
of these cases £ could be much larger than unity. For 
example, C ~ 990 if a = 11/6, M max = 0.1M e , and 
M m [ n = 10 M®. This value of M m ; n corresponds to 
the size of the smallest object that can survive radia- 
tion driven migration to the current age of the pulsar, 
as described in the analysis in Section 14.21 The value of 
-Mmax corresponds to the maximum size of an object that 
would remain undetectable in spectral analysis (see Sec- 
tion EXT) . 

Using cttoa = 45 M s for PSR B1937+21, Equation © 
constrains disk parameters according to 
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1 100 
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showing the clear trade off between the total mass in 
asteroids, the number of objects, and their distribution 
in semi-major axis and mass. For any configuration, the 
range of semi-major axes of the largest bodies can be con- 
strained by matching the spectral content of the residuals 
to simulations. 



relevant to PSR B1257+120 In this scenario a lower 
mass companion star evolves, over-flows its Roche lobe, 
and donates material to the neutron star, spinning it up 
from a radio-quiet object to an MSP with pulsed radio 
emission and a large particle wind. 

After (or perhaps during) recycling, the particle 
wind can ablate its com panion, as is the cas e in the 
PSR B 1957+20 system (|Fruchter et all [l98l and a 
growing number of similar pulsars. A number of the 
MSPs recently discovered in radio searches of Fermi 
;amma-ray po int sources are "black widow" objects 
Roberts! [201lT ) in which the MSP is ablating its for- 
mer mass-donor star. Additionally a growing number of 
MSPs are observed with low mass <C 0.1M© WD com- 
panions, with the MSP PSR J1719-1438 having nearly 
ablated its WD companion, transforming a star that was 
initially a large fraction of a solar mass in to an object 
with a mass of 10 -3 M Q (jBailes et al.ll2011h . 

The disk formation commences when the MSP disrupts 
the companion. While most of the material from the 
disrupted companion is accreted from the disk into the 
pulsar, conservation of angular momentum requires dis- 
persal of a fraction of the material over a wide region 
around the pulsar. The disk subsequently expands and 
cools enough to form the dust particles that seed the 
formation of larger objects. At this early time the disk 
may contain > 10 M s of dusty material extended over 
many AU (i.e. much more than in the present-day sys- 
tem, as we show below), well beyond the tidal disruption 
radius of the pulsar. Hi gh angular momentum systems 
may extend to 20 AU (jCurrie fe Hansenll2007ft . 

The high luminosity of PSR B1937+21 would pre- 
vent the formation of planetesimals and the survival 
of asteroids in the inner region of any disk. To esti- 
mate the radius where this happens, we assume ther- 
mal equilibrium of the material with input flux F- ln , 
and that the disk is optically thin. The input flux is 
dominated by the pulsar spin down luminosity L — 
Ak Ivv 25O2v0 where I is the neutron star moment 
of inertia and is assumed to be 10 45 g cm 3 . The pul- 
sar spin down luminosity comprises low-frequency Poynt- 
ing flux, a wind of relativistic particles, and high en- 
ergy magneotspheric radiation. The pulsar also emits 
thermal radiation that in the case of PSR B1937+21 
is negligible compared to this spin down luminosity 
(Wmai = 10- 4 - 7 L Q (i? NS /10 km) 2 (T/10 5 K) 4 , where 
i?NS is the neutron star radius and T is its temperature) . 

The flux incident on an object at orbital radius r is 



47rr 2 ' 



(11) 



3. ASTEROID FORMATION 

Several mechanisms have been proposed for the for- 
mation of planetary systems around pulsars. They 
include: (1) formation along with the p rogenitor 
(|Phinnev fe Hansen! H99l iPodsiadlowskil [19931] . (2) for- 
matio n from supernova fall back material ( Lin et aTl 
Il991h : or (3) form ation from material left over from 
the recycling phase (iNakamura fe Piranl ll991). The sec- 
ond is probably relevant to the disk around the magne- 
tar 4U 0142+61, and may also play a role in the vari- 
ability of the radio emiss ion observed in many pulsars 
(jCordes fe Shannonll200l . The third case is probably 



where e is the fraction of the total spin down loss rate 
that can be absorbed and gb is a factor that accounts 
for beaming of the pulsar's energy loss. The heating 
of the material depends on the fraction of the luminos- 
ity beamed toward the disk. In many models of pul- 
sar emission, the spin-down energy is beamed out of 
the open field line region around the magnetic poles 
(Goldrcic h fc Julian! I1969T) . However, recent numerical 
simulatio ns have suggested that the beaming is closer to 
isotropic (Spitkovskyl l201l] ). For the rest of the paper 

10 For an alternative explanation, sec Miller & Hamilton (2001). 
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present scaling laws based on the fiducial values e = 1 
and gb — 1. 

Assuming an asteroid radiates like a blackbody, the 
radiated flux is 



(12) 



where ctsb is the Stcfan-Boltzmann constant. If the as- 
teroid is in thermal equilibrium (i.e., ttR^F^ — 4nR 2 F out 
if the asteroid is spherical) , the equilibrium temperature 
T is 



T{r) = 1100 K (g b e 



,1/4 



L 



250L f . 



1/4 



/ r \-i/2 

(au) < 13 > 



suggesting that rocky material (melting point T mo i t ~ 
1100 K) can survive at r > 1 AU and a more metal-rich 
object (melting point » 1500 K) could survive to about 
1/3 of this distance. 

The formation of rocky macroscopic objects follows 
a path similar to that of planet formation in main 
sequence stellar system s and planet form ation in the 
PSR B1257+12 system (jHansen et alJl2009f ). where pro- 
toplanets undergo runaway growth until they reach a 
mass limited by their ability to gravitationally attract 
material from the surrounding surface mass density 
(|Lissauerlll987l ). 

Slower but inexorable growth of objects to planet sizes 
follows if perturbations induce orbit crossings. The rate 
of growth in this secondary phase depends on protoplanet 
masses, velocity dispersion, and the presence of gas to in- 
duce orbit crossings and accrete onto more massive plan- 
ets. 

Our proposal is that the evolution of the 
PSR B1937+21 differs from other known systems 
because its progenitor disk was too tenuous to form 
Earth-mass or larger planets. 

The largest isolation mass (the Hill mass) mu to which 
a protoplanet may grow during the runaway growth 
phase is the total mass co ntained within a n orbital radial 
range of several Hill radii (|Lissauerlll987l Equation 11): 



(4nBr 2 Z) 3 / 2 
(3M NS )V2 



10" 



^(r^So) 372 , (14) 



where B « 3.5 is a numerical constant and rpjj is the 
semi-major axis in AU, E is the surface mass density, 
and E is the surface mass density in g cm -2 . Without 
any means for inducing orbit crossing, objects do not 
feed well beyond their Hill spheres, resulting in a disk 
comprised of asteroid-mass objects. 

Planetesimals or protoplanets formed in such a disk 
would be too small to force gravitational interaction af- 
ter an initial runaway growth phase, and the disk would 
comprise many low-mass bodies in marginally stable or- 
bits. The disk would extend over a wider range of orbital 
radii than the solar system's asteroid belt because there 
are no planet-mass bodies to impose radial confinement. 

We suggest that disks of planetesimals may be a generic 
outcome of accretion driven spin-up and may potentially 
provide a noise floor for timing precision. The level of 
this floor would depend on the disk configuration, which 
itself evolves over the Gyr lifetime of the MSP. 

The surface density at the time of planetesimal forma- 
tion may have been larger if volatiles were present, and 



the Hill mass may be correlated with orbital radius and 
depend on the radial variation of E. Even in this case, 
our basic argument is unchanged. 

4. DISK STABILITY 

In this section, we investigate the stability of disks of 
asteroids to gravitational and radiative migration and 
show that these disks can persist to the present age of the 
pulsar. The spindown age of PSR B1937+21 is v/2\v\ « 
250 Myr and is an upper bound o n the elapsed time 
since the cessation of a ccretion (e.g. [Arzoumanian et all 
[19991 iWex etld][2000l) . iTaurisI (|2012j) includes a discus- 
sion on pre-accretion time scales and also suggests that 
spin down age likely overestimates the post-accretion age 
the system. Even in this latter case, the age is much 
longer than the expected runaway growth phase, but is 
short enough that the asteroidal disk we associate with 
timing residuals may still be evolving: slow coalescence 
may occur while the wind evaporat es the protoplanets 
that migrate to orbits near the NS (|Cordes fc Shannon] 
[20P.8t) . While direct collisions between disk objects are 
rare, gravitational stirring induced by a mass dispersion 
in the objects introduces a viscosity that will affect the 
lifetime of a massive system. Pulsar radiation clears the 
disk of its smallest objects. 

4.1. Gravitational Migration 

The arguments present ed here follow those presented in 
IHeng fc Tr emaind (|2010l ) , which examines the stability of 
planetesimals in a dynamically cold disk in which orbits 
do not cross. 

We first assume that all objects have a characteristic 
mass m, located in a disk of surface density E. In this 
case, the typical separation between objects is 



Ar 



m 
27rEr' 



(15) 



iChambers et al.l (|1996f ) examined the stability of plan- 
etary systems of up to 20 objects. They found that 
systems containing 20 objects became unstable after 
w 10 s Myr if the orbital separation satisfies 



Ar w 10r H . 



where 



r H = r 



\ 1/3 



3M NS 



(16) 



(17) 



is the Hill radius. 

IChambers et al.1 {l996) also found that the number of 
planets and the mass dispersion of the planets does not 
affect the stability time of the systems and argue that 
systems with orbital separations larger than 10 ru are 
likely stable because the time scale for interactions in- 
creases e xponentially with distan ce. We assume that the 
results of lChambers et a l. (1996) can be extended to sys- 
tems with a larger numbers of objects, and that objects 
spaced by 10 r# can persist for at least 250 Myr. 

For a belt of objects of mass m located between r\ and 
r2 the number of objects is therefore limited to 
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Substituting m 
we find 



N a < 



1000 



M hei t/N a and solving for N a , 
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We calculate the maximum-mass asteroid belt by us- 
ing the maximum allowed number of asteroids combined 

1/2 

with Equation (J7J), which implies Mb e it °c ctoa-^o • 
Numerically we obtain 
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4.2. Radiation Driven Migration 

Radiation forces will cause small objects to migrate 
through the system on time scales shorter than the age 
of PSR B1937+21. Here, we examine the effects of 
two different types of forces: radiation drag (Poynting- 
Robertson drag and an analogous effect due to the rela- 
tivistic pulsar wind) and Yarkovsky migration. 

The Poynting-Robertson effect causes small particles 
to migrate from the Earth into the Sun on a time 
scale t P R = 'jirpc 2 r 2 R/3L 672 yr pi£_ 4 r AU Z©/,L 
(|Burns et al.l (19 79), where p is the particle density in 
g cm~ 3 , i?_4 is the grain size in 10 -4 cm, tau is the or- 
bital radius in AU, and L is the luminosity. For a 10 /im 
grain orbiting PSR B1937+21 at 1 AU, the time for the 
grain to migrate to the pulsar is 



tp R w 2.7 yr(g b e) 



x (su)' 



lg cm 
L 



-3 



7? 



10" 



250L 



(21) 



Grains migrate very quickly in any compact disk sur- 
rounding the pulsar and little dust is expected close to 
the pulsar, unless there is a large reservoir of material 
further away. 

For much larger objects, the Yarkovsky effect will cause 
migration over the lifetime of the pulsar system. The di- 
urnal Yarkovsky effect is the net force on an illuminated 
object whose afternoon surface is hotter than its night- 
time surface. It can lead to an increase or decrease in 
semi-major axis depending on the lag angle of the hottest 
surface. In contrast, the seasonal Yarkovsky effect ap- 
plies to objects with spin axes tilted from their orbital 
plane, causing a net force along the spin axis that always 
decreases the semimajor orbital radius. 

Here we consider only the seasonal Yarkovsky effect, 
which yields a radial migration rate 



r«-5AU Gyr" 1 

r \-V2 / R a 



x 5b e 



1AU 



1 km 



L 



250 L, 



, (22) 



where we have re- sc aled Equation (A4) of 
ICordes fc Shannon! ((2008) using fiducial values bet- 



ter matched to the disk configuration under analysis 
here. 

We define a critical asteroid size i? yar such that all 
objects smaller than i? yar migrate inward from a distance 
r to the pulsar magnetosphere(3 in a time T. Scaling to 
PSR B1937+21, we find, 

i? yar = 5.1 km 9b e ^ (JLS) {^y^ 

The Yarkovsky effect therefore can cause objects up 
to ~ 5 km to migrate inward over the Gyr lifetime of 
many MSPs and presumably be evaporated well outside 
the light-cylinder radius. For spherical asteroids with 
densities of 1 g cm -3 , i? yar corresponds to objects with 
masses of 10 -10 M®. The Yarkovsky effect cannot 
cause larger objects than this, or much more distant than 
1 AU, to migrate appreciably over the lifetime of the 
system. 

We conclude that a disk extending beyond ss 1 AU, 
populated by objects with sizes larger than ~ 5 km 
would survive radiative effects up to the present day for 
PSR B1937+21 and thus could induce timing variations 
at the present epoch. 

5. PLAUSIBLE DISK CONFIGURATIONS 

The residual TOAs for PSR B1937+21 are domi- 
nated by low-frequency components that are nevertheless 
spread over a fairly broad band. To find asteroid systems 
that reproduce the frequency content of the data, sim- 
ulations of a wide range of asteroid belt configurations 
were conducted. Each asteroid was randomly assigned 
an orbital radius according to a probability distribution 
that is flat in radius (hence decreasing in asteroid sur- 
face density with radius) between an inner radius r% and 
an outer radius ra. Within its orbit, each asteroid was 
given a random initial phase. The masses of the objects 
were selected from a power law mass distribution. The 
asteroids were modeled to be non-interacting, which is 
an appropriate approximation because the evolution of 
the orbits over the 26 yr observing span (corresponding 
to at most a few orbits) is small. Additionally, all ob- 
jects within the simulated disks were required to satisfy 
the dynamical and thermal migration criteria defined in 
Equations ([l~6"l) and (|23|) . respectively. 

Random placement of objects alone was found to 
be computationally inefficient at producing dynamically 
stable disks when the number of objects was greater 
than « 500, because the probability of an object be- 
ing placed within 10 rjj of another object is large. To 
increase computational efficiency, we edited out objects 
that were in dynamically unstable orbits. Disks were it- 
eratively populated disks with an increasing numbers of 
objects. If any object was initially placed in an orbit 
within 10 Hill radii of another object it would be deleted 
from the disk. Consequently, simulated disks contained 
fewer objects than initially specified. We found that it 
was also necessary to edit simulations that contained de- 
tectable individual sources. The minimum detectable as- 

11 Objects that have migrated inward will be evaporated well 
before they reach the magnetosphere of the MSP (defined as the 
light cylinder radius, tlc = cP/2n S3 74 km.) and thus are not ex- 
pecte d to produce any intermi ttency effects or torque fluctuations 
(e.g., [Cordcs & Shannon^ 
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teroid mass is discussed in Section [7. 1.1 1 Based on anal- 
ysis of that section, we restricted objects to have masses 
m < 5 x lCT 4 Af e (77l yr) 2 for / > 1/26 yr" 1 . 

At each epoch (corresponding to the observation times 
in the actual data), the cumulative orbital effect from 
all asteroids was calculated. To this, a random white- 
noise contribution of rms amplitude equal to the formal 
TOA uncertainty was added. The residual TOAs were 
found by subtracting from the simulated TOAs a least- 
squares fit comprising quadratic terms that model pulsar 
spin down. To verify that the simulation code worked, 
we generated systems comprising one or a few objects 
with known masses. We found that the amplitude of 
the TOA perturbations agreed with the expressions in 
Equation ((TJ). 

The real and simulated data were compared using max- 
imum entropy power spectra of the timing residuals. The 
quality of the match was assessed using a cost function, 



C = ^[ln(4°)/5 



(24) 



where and are the power spectra of the observed 
and simulated data sets, respectively. We searched for 
the minimum-cost combinations of parameters by using 
a brute-force search over a grid of values for parameters 
that are displayed in Table [1] For each unique combi- 
nation of parameters, 100 stable asteroid belt configu- 
rations were simulated. We only simulated disks com- 
prising N a < 4000 objects with a minimum mass much 
larger than the minimum mass object that can survive. 
For power law mass distributions with a < 2, the mass 
dispersion (m 2 ), and hence cttoa (see Equation 0]) are 
dominated by the largest objects, and (for our purposes) 
it is not necessary to simulate the smallest objects. This 
was done primarily to reduce the computational burden. 

No unique combination of parameters was found that 
minimized the cost function. This is highlighted by two 
of the best-matching asteroid belt configurations, which 
are displayed in Figure |H One corresponds to an asteroid 
belt with N a = 80 objects, total mastQ M = 0.04M e , 
and a relatively flat mass distribution (a = 1.1), while 
the second has 1609 objects, a total mass of M = 0.06M® 
with a steeper mass distribution (a — 1.5). Both belts 
extend to an orbital radius of approximately 15 AU and 
therefore contain objects orbiting at periods longer than 
the data span length. 

The number of asteroids can vary in the best-matching 
configurations from a few tens to a few thousands, but 
the total asteroidal mass is found to approximately sat- 
isfy Equation (fl0|) . 



6. COMPARISON WITH OTHER PLANETARY SYSTEMS 

iCurrie fc Hansen! (|2007D and lHansen et all ([20091 ) 
modeled the formation of solid objects in disks formed 
by either supernova fallback or tidal disruptions of com- 
panion white dwarf stars with an emphasis on modeling 

12 The total disk masses do not match the input grid-search 
values listed in Table [T] for two reasons. Firstly, asteroids were 
randomly assigned masses, so agreement with the injected value is 
only expected in the ensemble average. Variation about this mean 
is expected, especially in cases where the mass distribution has 
a steep negative spectral exponent. Secondly, as justified above, 
dynamically unstable objects were culled from systems, reducing 
disk masses. 



the PSR B1257+12 planetary system. The predominant 
difference between the two types of disks is their initial 
angular momentum, with the tidal disruption disk con- 
taining larger initial angular momentum. They found 
that both kinds of disks produced planets in many-AU 
orbits, much larger than the 0.2 to 0.5 AU orbits of the 
planets around PSR B1257+12. The tidal disruption 
planets were found to extend to orbital radii > 20 AU. 
These disks compare favorably to our best-fit models 
which have objects orbiting at distances of ~ 20 AU. 

Comparison with our maximal-mass solution implies 
that at 1 AU, the present surface density is S w 
0.013 g cm" 2 , smaller by three orders of magnitude 
than the minimum- mass surface density (i.e., excluding 
volatiles) thought to be pr esent in the solar system a t the 
end of the runaway phase (jLissauer fc Stewartil993[ ) , and 
smaller by > 10 than in the planet formin g region near 
the planet pulsar system PSR B1257+12 feudenl 119931: 
ICurrie fc Hansen! 120071 : lHansen et all |2009J). However, 
our estimate of the surface density may underestimate 
the initial surface density due to the removal of material 
caused by the effects discussed in Section S) 

7. TESTS OF THE ASTEROID BELT MODEL 

A drawback of an asteroid belt cause for the RN ob- 
served in the timing residuals of PSR B1937+21 is that 
it is difficult, though not impossible, to test. In the fol- 
lowing subsections we first highlight tests that can be 
conducted using high precision pulsar timing, and then 
discuss tests generic to detecting circumpulsar objects. 

7.1. Pulsar Timing Tests 
7.1.1. Detecting the Periodicity of Individual Objects 

The most direct way to confirm the asteroid belt hy- 
pothesis is to detect one or more individual objects by 
identifying features in the power spectrum that are sta- 
tistically significant. The power spectral density formed 
from the TV-point discrete Fourier transform Xk of dis- 
cretely sampled data is 



T ~ 



(25) 



where T is the observing span. 

Including additive white noise with the timing pertur- 
bations from a collection of N a non-interacting asteroids, 
the power spectrum has an ensemble average 

(Pk) = (Pk,a) + (Pk, W ), (26) 

where (Pk,w) is the mean white-noise spectrum, 

Tal 



Pk,w) = 



N ' 



(27) 



where a w is the rms white noise error. 
The mean asteroid-induced spectrum is 



(P 



T 



N„ 



2 ^ 



a\ [sincd 2 7r(fc/7V- f a At) 



sincd 2 7r(fc/A + / Q At)] 



(28) 



where we have used the discrete "sine" function defined 
as sincd(a;) = sin Nx/N smx, which is unity for x = 
0. In the null hypothesis where there is no signal, the 
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Table 1 

Grid search parameters 



Parameter 



Description 



Values 



M a Total mass (Af®) 

N a Number of asteroids 

a min Inner radius of asteroid belt (AU) 

imax Outer radius of asteroid belt (AU) 

a Spectral index of mass distribution 

A^min Minimum mass of asteroid (Mj ) 



0.001, 0.003, 0.01, 0.03, 0.1, 0.3 

3, 10, 30, 100, 300, 1000, 2000, 4000 

0.5, 1.0, 3.0, 10.0, 30.0 

1.0, 3.0, 10.0, 30.0. 100.0 

0.1, 0.5, 0.9. 1.1, 1.25, 1.5, 1.67, 1.83 

lO" 8 



Note. — Grid Search Parameters. For each unique set of parameters, 100 stable asteroid 
belt configurations were simulat ed and compared to the observed residuals using the cost 
function expressed in Equation (l-'4[ i. 
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Figure 4. Best-fit asteroid belts and their signatures in residual TOAs. In the top row (panels a-c) we show a configuration comprising 
80 objects with a total mass of 0.04 M®. In the bottom row (panels d-f) we show a configuration comprising 1609 asteroids with a total 
mass of 0.06 Mq. Leftmost panels (a and d): Asteroid belt configuration at initial time in the simulation. Each asteroid is represented 
by a filled circle, with the size of the circle proportional to the radius of the asteroid. Open circles highlight objects with small masses. 
Center panels (b and e): Residual TOAs associated with the asteroid belt. Rightmost panels (c and f ): Maximum entropy power spectra 
for the simulated residuals (thick lines) and the observed time series (thin lines). 



power spectral estimates have a x 2 distribution with two 
degrees of freedom. Defining a normalized spectrum to 
be 

Pk 



{Pk, 



(29) 



for the null case, the probability that any single value Sk 
exceeds a value S is 

P{S k > S} = exp(-S). (30) 

When M = N/2 ^> 1 spectral amplitudes are tested and 
we require the mean number of false positives to be less 
than one, 



we obtain 



Mexp(-S) < 1, 
S > InM. 



(31) 
(32) 



More generally, the probability that one or more out of 
M spectral values exceed S is the false-alarm probability 

Pm 



Pm = P>i(> S) = 1 - [1 - exp(-S) 



(33) 



An asteroid signal would need to be larger than S for 
it to not be considered a false-positive; equivalently, we 
can say that the null hypothesis can be rejected with 
probability 1 — pm ■ This implies that 



S > 

or, for p M < 1, 



In 1-(1-pm) 



S > In M/pm- 



l/M 



(34) 



(35) 



Detecting Periodicities in a Sparse Asteroid Belt — We first 
consider a sparse asteroid disk that contains only a few 



Circumpulsar Asteroids 



9 



objects that produce spectral lines in Equation (j2"5)l that 
do not overlap and thus potentially can be identified. In 
this case, detectability is limited solely by the additive 
white noise. For an asteroid signal that maximizes at k a , 
this requirement becomes approximately 



P ka > (Pk,os)HM/p M )- 



(36) 



Substituting from Equation (|28[) and using -/\ ff 
Ta^/N the asteroid signal must satisfy 

1/2 



02 > cr w \ — 



1 ( M \ 
In 

Pm 



1 



1/2 



(37) 



where a?, is defined in Equation [TJ 

The implied minimum detectable asteroid mass is then 



V2cM NS 



asiniF(P olh ,T) 
M 



x 



In- 



Pm 




(38) 



Substituting values appropriate for PSR B1937+21 in- 
cluding a neutron star mass of 1.4 Mq, using a false- 
alarm probability pu — 1/M (i.e. a threshold corre- 
sponding to a single false alarm), and M = N/2, the 
minimum detectable asteroid mass is 



10- 49 M ( 



l» n 



smiF{P orhl T) 



V100 ns/ V N 



1/2 



Porb 



-2/3 



2 In TV — 1 
21nl0 3 - 1 



1/2 



(39) 



Individual objects of approximately one third the mass 
of the main belt asteroid Ceres can be detected in the 
Fourier transform if they are orbiting the pulsar at 1 AU, 
and only white noise is present in the observations. 

Resolving spectral lines associated with a large num- 
ber of objects is more difficult, however. If the aster- 
oids comprise a marginally stable and maximally packed 
belt with separations Ar/r rj 10(m/3MNs) 1 ^ 3 , the fre- 
quency resolution Sf = T _1 must be smaller than the 
separation of spectral lines produced by neighboring as- 
teroids Sf/f = (3/2) Ar/r, which, for a 1.4M neutron 
star works out to 



T > 161 A0-M fl 



P 



orb 



1/3 



(40) 



For data spans of 10 to 20 years, the frequency resolution 
will be sufficient to resolve the spectral line from an in- 
dividual asteroid only if the disk is much less populated 
than a maximally packed disk or if the asteroid has a 
mass much larger than the average. 

Detecting Periodicities in a Dense Asteroid Belt — Even if 
the cumulative effect of the asteroid belt produces a red 
noise-like signature, it may still be possible to isolate 
the periodicity from the largest object or objects in the 
disk. Using the PSD in Figure [5] as a basis, we model the 
total noise background from which an individual asteroid 
needs to be distinguished by a combination of white noise 
and a steep red component: 



P k ,os{k = fT) = 
0.015 /is 2 yr 



lyr 



-5 



0.015 ^s 2 yr. (41) 



The minimum detectable mass can be found by substi- 
tuting Equations (|4"Tj) and (|2"5)) into Equation (|3"6f and 
solving for m. In Figure [5] we show a plot of the mini- 
mum detectable mass as a function of orbital period P or b- 
When P or b ^ 1 yr the sensitivity is poor because of large 
levels of red noise. When P or b -C 1 yr, white noise dom- 
inates and sensitivity is poor because the reflex motion 
per unit mass is diminished. Over a wide range of or- 
bital periods we can rule out any single object with a 
mass < 0.05M®. In Figure [5] we also show the distribu- 
tion of objects in mass and frequency from one of our 
best simulations. As expected, none of the objects ex- 
ceed the threshold for detectability in a search for their 
periodicity. 




orb 



(yr 



Figure 5. Minimum detectable mass of an asteroidal compan- 
ion to PSR B1937+21 for a dense asteroid belt. The solid lines 
show the sensitivity assuming tha t th e residuals contain red noise 
of the form presented in Equation 1411 1 , which is a model consistent 
with the observed red noise. From bottom to top, the lines repre- 
sent false alarm probabilities of p = 0.33, 0.05, and 10 -4,2 , which 
roughly correspond to 1 — <r, 2 — a, and 5 — cr limits. The dashed 
curves show the sensitivity assuming that the only source of noise 
is white noise for the same set of false alarm probabilities. The 
filled circles denote the three planets orbiting th e millisecond pul- 
sar PSR B1257+12 (IKonacki fc W olszczan 2003), which all would 
have been detected if present around PSR B1937+21. The open 
circles denote show the masses of the objects the best-fitting config- 
urations displayed in panels d — f of Figure [4] The vertical dashed 
line shows the data span length. The change in slope and loss of 
sensitivity for P or b > 20 yr is associated with the polynomial fit 
for pulsar spin down. The narrow spikes at 0.5 and 1 yr correspond 
to the fit for parallax and proper motion, respectively. 



7.1.2. Stationarity in Time Series 

If individual asteroids cannot be discriminated, 
the statistics of their collective behavior can be 
used to distinguish the asteroid model from intrin- 
sic spin noise. Spin noise is manifested as non- 
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stationarit-JTCl in the res iduals (jCordes fc Do wns 1985; 
IShannon fc Cordesl [20loh that can be modeled as ran- 
dom walks in spin phase, spin frequency and fre- 
quency derivative, and other similar processes hav- 
ing power spectra that are power-law in form. For 
these processes respectively, the rms residuals scale 
as cr R N, 2 (T) oc T 7 , with 7 = 1/2,3/2,5/2. Most 
CPs show consistency with the frequency and fre- 
quency derivative spin noise, a superposition of the 
two, or rando m walks combined with occasional dis- 
crete jumps (jCordes fc Downslll985t ID'Alessandro et al.l 
I1995D . More recently IShannon fc Cordesl (|2010D found 
that for rotation-powered pulsars CTRN,2(r) oc T 2 - 0±0 - 2 . 
The dependence on the data span T is a clear signature 
of non-stationary statistics. 

In the asteroid model, for time scales shorter than the 
largest orbital period (i.e, T < P rb,max), crn,2 increases 
with T because larger fractions of the perturbations asso- 
ciated with the longest period objects are manifested in 
the post fit residuals. However, long term stationarity is 
expected from asteroid belts for T 3> -Porb, max- The tran- 
sition from non-stationary to stationary statistics can be 
explored by applying a variety of diagnostics to the resid- 
ual TOAs, including structure functions, spectral estima- 
tors, and <trn,2(T') itself. If non-stationarity persists for 
long periods of time (much greater than 25 years) , either 
an asteroid belt model would need to be developed that 
includes objects in very wide, stable orbits, or the as- 
teroid belt origin for observed TOA variations could be 
ruled out. 

7.1.3. No Asteroids Inside the Rock Line 

As mentioned earlier, the spin-down luminosity of 
PSR B1937+21 is large enough to heat asteroids to the 
large temperatures needed to evaporate them (c.f., Equa- 
tion [131). For PSR B1937+21, asteroids composed of 
rocky material would not be expected inside 0.5 AU 
(P orb « 0.4 yr), assuming a melting temperature of 
w 1500 K and complete absorption of the pulsar spin 
down luminosity (e = 1) and isotropic beaming (g^ = 1). 
This rock line is analogous to the water/ice line discussed 
in planetary and exoplanetary science. Asteroids com- 
posed of ice and other refractory materials would not be 
expected inside a radius of 16 AU (P or b ~ 14 yr), as- 
suming a melting temperature of 300 K and complete 
absorption of isotropically beamed pulsar spin down lu- 
minosity. The absence of asteroids inside the rock line 
would be manifested as a paucity of spectral power at 
fluctuation frequencies / > l/-P rb,rock> where P rb,rock is 
the orbital period at the rock line. The spectra in Fig- 
ures Q]and[4]are consistent with this picture. In order to 
detect the presence of red noise in this region, it would be 
necessary to reduce radiometer noise, by observing with 
a higher gain system, a wider bandwidth, or a longer 
integration time. 



13 We note for some canonical pulsars, it appears that there is a 
sub-dominant quasipe riodic sourc e of RN attributed to magneto- 
spheric activity (Lvnc ct al. 2010). Wc refer the reader to Section 
17. 1.61 for further discussion. 



7.1.4. Differences in Red Noise between Isolated and Binary 

MSPs 

Relative to binary MSPs, the processes that form iso- 
lated MSPs may provide initially more massive and ex- 
tended disks of material. This is plausible because iso- 
lated MSPs have presumably completely disrupted their 
companion and polluted their environment with a larger 
mass of material with a wid er range of angular momen- 
tum (jCurrie fc Hansen 2007). In this scenario, as a pop- 
ulation, isolated MSPs would possess larger and more ra- 
dially extended disks, which would be expected to have 
larger objects in wider orbits and would have larger levels 
of red noise, as expressed in Equation (J7J). 

It should be noted that binary systems themselves 
can potentially possess stable circumbinary disks because 
only the inner regions of the disk are strongly torqued 
by the binary pair. Indeed, long lived debris disks 
have been observed arou nd main sequence binary sys- 
tems (|TriHing et al.l [2007h ■ In numerical simulat i ons o f 
late-stage planet formation. iQuintana fc Lissauerl (2006) 
found that planet formation surrounding compact binary 
star systems was statistically identical to planet forma- 
tion around isolated stars, provided the binary systems 
were c ompact and in circular orbits. iHolman fc Wiegertl 
(1999) examined the stability of single planet orbits 
around binary systems, in which the binary has a wide 
range of eccentricities and mass ratios. Through numeri- 
cal simulation, they derived an empirical relationship for 
the innermost stable orbit as a function of semi-major 
axis, eccentricity, and mass ratio of the binary. For mass 
ratios between 1:1 and 9:1, and systems with eccentrici- 
ties between and 0.7, they found that the planets were 
stable with orbits that were a factor of four greater than 
the semi-major axis of the binary. However, simulations 
were conducted for only 10 orbits, which is more than 
5 orders of magnitude shorter than the time required for 
MSP systems. Because NS-WD binary orbits are typi- 
cally very compact (P or b few hr to 80 d), asteroids or- 
biting at > 1 AU are not strongly affected by the inner 
binary pair. 

Disks in principle could persist around double neutron 
star (NS-NS) systems as well. However, pulsars in NS-NS 
binaries do not exhibit the timing precision that MSPs, 
even without asteroidal effects. Firstly, timing residuals 
from NS-NS binaries are larger because pulse widths and 
there are fewer pulses per unit time interval, so pulse- 
phase jitter is larger. Secondly there appears to be timing 
noise and intrinsic spin instabilities in some of the objects 
(|Weisberg et al.l [20101) . Therefore the perturbations of 
an asteroid belt that is allowed around B1937+21 would 
be masked by the other timing effects. An asteroid belt 
would have to be accordingly more massive or have a 
small number of more massive individual asteroids for 
there to be a competing effect from asteroids. 

At present, there is sparse and perhaps only cir- 
cumstantial evidence that MSPs with massive compan- 
ions are intrinsically more spin stable than isolated 
MSPs. Indeed, neither the PSR B1257+12 system 
nor PSR B1937+21 have stellar companions. By con- 
trast, the MSPs with the best timing stability have 
WD compan ions (PSRs J0437-4715, J1713+0747, an d 
J1909-3744: lVerbiest et"al|[200l iDemorest et al.ll20ll . 

However there are a number of isolated MSPs that do 
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not show the cubic signature of strong RN at the cur- 
rent levels of timin g precision (e.g., PSR J2124— 3358, 
iVerbiest et al.l [20090 . Evidence for RN in more isolated 
MSPs would be necessary to show a definitive difference. 
Even then, careful analysis would be required to conclu- 
sively identify that the correlation is due to an asteroid- 
belt origin and not other plausible origins for red noise. 

7.1.5. Asteroid Belt Size and Composition Correlated with 

MSP Age 

Asteroid belts may evolve slowly over the age of the 
pulsar. Radiation, gravitational migration, and occa- 
sional collisions will cause objects to migrate inward to- 
ward the pulsar, particularly if the disk is more dense 
than minimally required. As material migrates inward, 
the reflex motion of the NS is diminished. Disks would 
be expected to be more compact, and (potentially much) 
smaller in mass around the oldest MSPs. 

7.1.6. Residuals Uncorrelated with Magnetosphenc Activity 

iLvne et al.l (|2010j ) found that some of the RN in some 
CPs is associated with changes in state between dis- 
crete values of the spin down rate, and that the 
states are associated with different average profile shapes. 
While there is no evidence for this process occurring in 
PSR B1937+21, if all of the RN was associated with this 
type of behavior, an asteroid belt origin would be ruled 
out, because reflex motion would neither cause transi- 
tions between discrete states of v nor cause the temporal 
variability of pulse profiles. As discussed earlier, the level 
of RN f or this pulsar is consistent with what is observed 

m CPs dsE 

annon fc Cordesl[20T0l) . However, given pul- 
sar to pulsar variability in the levels of RN in CPs, some 
fraction of the observed RN in this PSR B1937+21 may 
still be associated with reflex motion. 

7.1.7. Variations in Dispersion Measure 

Another possibility is to search for variations in dis- 
persion measure (DM) due to evaporation of asteroidal 
material by the pulsar's relativistic wind. Whether DM 
variations are detectable depends on the evaporation rate 
and on how fast thermal gas is entrained in the relativis- 
tic flow. Relativistic gas will produce a smaller timing 
perturbation than thermal gas with the same density. 
Some fraction of the total measured DM from B1937+21 
could be associated with thermal gas near the pulsar, 
most likely just outside the bow shock that, for this 
object would be 0.04 p c from the pulsar using ap pro- 
priate scaling laws (e.g. iChatteriee fe Cordesl 120021 ) for 
an ambient density of 1 H atom cm -3 and a nominal 
three-dimensional space velocity of 100 km s -1 . How- 
ever, evaporation events much closer to the pulsar might 
produce short-lived increases in DM, but whether these 
occur or not is related to the location of the rock line, 
as discussed above, and whether small asteroids migrate 
to lower radii where they evaporate more quickly. Re- 
gardless, variations in DM could provide a handle on 
circumpulsar debris. Any timing variations would need 
to be distinguished from interstellar variations, which 
produce ~ 2 /us at 1 . 4 GH z over time scales of years 
(|Rama chandra n~et al.l 120061). and ~ 100 n s variations 
over time scales of 10 days ([You et al.ll2007ft . DM varia- 
tions from B1937+21 are consistent with an interstellar 



origin because they conform to the expectations for a 
Kolmogorov medium combined with discrete structures. 
Together these produce flux-density vari ations that are 
anti-c orrelated with timing variations (jLestrade et al.1 
1998). Flux density variations would not be expected 
from evaporation events near the pulsar that conform 
to the observed anti-correlation, which is caused by re- 
fraction in a cold plasma. While DM-induced timing 
variations are larger at lower frequencies, scaling propor- 
tional to v~ 2 , multi-path propagation delays associated 
with the ISM scale even more strongly (oc v~ A ) and they 
can limit the timing precision necess ary to determine 
the d ispersion measure at each epoch ((Foster fc Cordesl 
1990). Thus, while searching for DM variations from cir- 
cumpulsar material may be difficult for PSR B1937+21, 
other MSPs, especially those that are nearer the solar 
system and thus have smaller ISM contributions, may be 
fruitful targets for low -frequency arrays, such as LOFAR 
(jStappers et al.ll201lD . 

7.2. Other Tests Asteroid Belts around Pulsars 
7.2.1. Infrared Emission from Debris 

Any debris around neutron stars would be passively 
heated by the neutron star's thermal and non-thermal 
emission to temperatures of hundreds of Kelvin, radi- 
ating predominantly in the infrared. For debris to be 
detectable, the disk would be required to have a large ef- 
fective area. Infrared observations of other pulsars have 
placed limits of many Earth masses on circumpulsa r dust 
disks (jLazio fc Fisc her 20041: IBrvden eFall 12006). Be- 
cause of the large d istance to PSR B1937+21 (5^ kpc, 
IVerbiest et al.ll2~012l ). tens to hundreds of hours of tele- 
scope time with a JWST-like telescop e would be required 
to de tect an asteroid-belt like disk (jCordes fc Shannonl 
2008). The lack of a population of small radiating parti- 
cles with a cumulative large effective area (removed from 
the system due to Poynting- Robertson drag) further in- 
hibits the detection of the disk. In addition, definitive as- 
sociation of infrared emission would be severely hindered 
by confusion with coincident infrared sources because the 
pulsar is located in the Galactic plane. 

7.2.2. Reflection of Radio Emission off of Asteroids 

If the pulsar beam intersects the asteroid belt, radio 
emission will reflect off of any material with siz e greater 
than the wavelength of the radio emission ((Phillips 
119931 : iCordes fc Shannonl l2008h . Reflected radio emis- 
sion can be distinguished from magnetospheric emis- 
sion through polarimetric observations because reflected 
emission would be depolarized by scattering off of the 
rough asteroid surfaces, while the magnetospheric emis- 
sion would show substantial levels of polarization. De- 
tection of this emission is challenging with single-dish in- 
strumentation because the faint off-pulse emission cannot 
be distinguished from receiver noise and other sources 
within the comparatively larger beam. In the past, ra- 
dio interferometers have not been capable of producing 
the fast-dump imaging necessary to distinguish pulsed 
point source emission from non-pulsed reflection emis- 
sion. However, new instrumentation at the Australia 
Telescope Compact Array, Giant Metrewave Radio Tele- 
scope, Jansky Very Large Array, MeerKAT (under con- 
struction), and the planned Square Kilometre Array will 
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be capable of fast-dump imaging and it will be possible 
to identify off pulse emission. Recently, using the Giant 
Metrewave Radio Telesco pe, off-pulse emission h as been 
associated with two CPs (|Basu et al.ll201ll . l2012t ). How- 
ever, this emission has been attributed with the pulsars' 
magnetospheres . 



8. IMPLICATIONS: TIMING PRECISION AND 
GRAVITATIONAL WAVE DETECTION 

If MSP timing precision is improved further it will be 
possible to detect gravitational waves (GWs) by observ- 
ing time of arrival (TOA) variations that are correlated 
amon g a set of MSPs comprising a pulsar timing array 
(P^ jDetweTlel1[l^7§lF"oster fe Backerll990l:|Jen"et et all 
120051) . The most plausible source of GWs is a stochas- 
tic background of gravitational waves ( GWB) associ- 
ated with massive black hole binaries (e.g.. ISesana et all 
120081) . 

There are three ways to improve the sensitivity of a 
PTA to gravitational waves: observe a larger number 
of pulsars, make more precise estimates of TOAs, or 
observe over a longer timing baseline. As long as tim- 
ing error is dominated by the finite signal to noise ratio 
of the observations (i.e., radiometer noise), timing pre- 
cision can be improved by observing the pulsars more 
frequently, with longer integration times, wider band- 
widths, and higher gain radio telescopes. However there 
are ma ny reasons to expect limitations to these improve - 
ments (IQslowski et all 12011 IShannon fe Cordesl |20T1 . 
including the presence of asteroid belts around MSPs. 

In the previous sections, we have demonstrated that 
MSP systems that contain a population of large, rocky 
objects possess poorer timing residuals than those that 
do not. Assuming the disks consist of more than a few 
objects it is nearly impossible to mitigate the effects of 
these disks in residual TOAs. Most obviously, orbital 
noise cannot be reduced by simply performing more sen- 
sitive observations, e.g., by using more sensitive backend 
instrumentation or a higher gain telescope. Secondly and 
more alarmingly, the signature of orbital noise is similar 
to that imparted by the GWB on pulsar TOAs. As a 
result, for a given rms signal strength, the tolerance is 
much lower for orbital noise than white noise sources: 
reflex motion inducing 20 ns rms variations over 5 to 
10 yr wil l significantly reduce the s ensitivity of a PTA to 
a GWB (IShannon fe Cordesll2010l) . 

In the presence of red noise, there are two ways to 
improve PTA sensitivity to GWs. A modest increase 
in sensitivity can be achieved by increasing observing 
throughput, i.e, by observing some combination of higher 
cadence and longer per-epoch observations. The only 
way to significantly increase sensitivity is to incorporate 
additional MSPs with comparable (or superior) stability 
into PTA observations, such as pulsars with lower lev- 
els of intrinsic red noise or smaller (or absent) asteroid 
disks. Indeed there are a handful of pulsars that are sen- 
sitive to the GWB within its pred icted amplitude range 
and are absent of strong red noise ([Demorest et alj[2013t 
iManchester et all 12012;. For the purposes of detection, 
it is imperative to identify other systems that show sim- 
il ar stability. These implica tions are discussed in detail 
in lCordes fe Shannon! (f20ll . 



9. CONCLUSIONS 

We have shown that asteroid belts can form and per- 
sist around MSPs, and these systems can affect the tim- 
ing precision at the ns to /xs level. The residual RN 
in observations of PSR B1937-I-21 is consistent with the 
presence of an asteroid belt of mass ~ O.O5Af0. The 
mass and distribution of objects of objects in the disk 
is constrained by the timing variations. However, addi- 
tional constraints can be placed because the disk must be 
dynamically stable and asteroids need need to withstand 
the strong radiation field emitted by the pulsar. In Table 
[2] we list constraints on the disk configuration derived in 
this paper, including the minimum and maximum mass, 
the mass and spatial distribution of the asteroids, and 
the number of objects. 

The existence of planets around two pulsars, a debris 
disk around another, and a possible asteroid belt around 
PSR B1937+21 have several implications. Protoplane- 
tary disks around pulsars must cover a wide range in 
density and total mass in order to have planets form in 
some cases but not in others. Disks of lower mass may 
exist around other MSPs, but not have been detected 
because of insufficient timing precision. 

Although it is difficult to test the asteroid model, it 
may be equally difficult to argue that absolutely no de- 
bris would be left in the system once accretion driven 
spinup and evaporation of the companion were com- 
plete. Confirmation of the orbital noise interpretation for 
PSR B1937+21 may rely on finding planets and asteroids 
around other pulsars and demonstrating that they are a 
frequent outcome of processes that lead to recycled pul- 
sars that presently do not have a single large companion. 

If disks are populated by more than a few objects, 
it is difficult to identify individual objects and correct 
for their TOA perturbations. Unless PSR B1937+21 is 
unique, high precision timing observations of other MSPs 
may show unmitigatable, correlated orbital noise associ- 
ated with circumpulsar disks. More MSPs need to be 
investigated in detail to search for red noise signatures 
consistent with astroid disks. 
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Table 2 

Constraints on Disk Configuration 



Constraint 


Mechanism 


Reference 


Comment 


Mass Distribution (obs) 
Maximum Asteroid Mass (obs) 


Reflex Motion 
Reflex Motion 


Eqs. (4)-(10), (20) 
Figure 5 


Dependent on asteroid mass distribution (factor f in Eq. 8) 
Constrained by power spectral analysis to Af max ft 0.05 


Minimum Asteroid Size 
Minimum Orbital Radius 
Orbital Separation 
Maximum Number of Objects 


Yarkovsky Migration 
Radiative Heating 
Dynamical Stability 
Dynamical Stability 


Eq. (23) 
Eq. (13) 
Eqs. (16)-(17) 
Eqs.(18)-(20) 


R mirL ft 5 km, M min ft 10- lu Af e at 1 AU 
a min « 0.5 AU, P orb min ft 0.3 yr if T melt = 1500 K 
Ar > lOrjj (extrapolated from Chambers et al. 1996) 
Dependent on minimum orbital separation 



Note. — Constraints on asteroid belt configurations derived in this paper. For each constraint we list the mechanism that provides the 
constraint, the relevant equation or figures in the text, and a brief description of the magnitude of the constraint or assumptions required to 
set its magnitude. Constraints that are derived from analysis of the time series are listed as (obs). 



APPENDIX 

A. OBSERVATIONS OF B1937+21 

In this appendix, we outline the methods used to form the residual times series displayed in Figure [T] and the power 
spectrum in Figure [5] that are the basis for the analysis in this paper. Observations are from the Arecibo, Effelsberg, 
Nangay, and Westerbork telescopes. In Table [31 for each telescope we list the MJD range of observations, observing 
span T, nominal observing frequency ^, observing backend, the typical uncertainty in TOA, and the total number 
of TOAs ./Vtoa- In this table, we also give a reference that describes the instrumentation. All of the observations 
from the Effelsberg and Nangay telescopes were conducted in bands close to 1400 MHz. Observations at Arecibo were 
conducted in two widely separated frequency bands, while observations at Westerbork were conducted in three widely 
separated bands. All observations were corrected for known offsets in the observatory clocks, including a clock offset 
recently found in Westerbork data (G. Janssen et al., in preparation). The TEMP02 package (H obbs et al.ll2006l ) was 
used to analyze the TOAs. In the next two subsections, we describe how the timing solution used to for this analysis 
was made. 

A.l. Correcting for Propagation Effects 

It is well known that it is necessary to correct for interstellar propagation to achieve the best timing precision. 
Because the refractive index of the interstellar medium is frequency dependent, propagating radio waves at different 
frequencies arrive at different times and travel along dif ferent paths. This is m anifested as a number of perturbations 
of the observed signal relative to the emitted signal (e.g. lFoster fc Cordeslll99C)f h At minimum it is necessary to correct 
for the group delay associated with the total electron content of the line of sight, the latter referred to as the dispersion 
measure (DM) in the parlance of pulsar astronomy. 

DM variations i n the Arecibo dat a were corrected using the published DM time series that accompany the dual 
frequency TOAs (jKaspi et al.lll994D . For the earliest Nanqay o bservations (MJD < 52579), DM variations were 
corrected using published DM values (jRamachandran et alJ l2006f ) that were partially derived from the observations 
presented here. For later observ ations, the multi-frequ ency Westerbork observations were used to model DM variations 
using a ninth-order polynomial (jLazaridis et al.ir2011[ ). While the DM variations are significant, their contribution to 
the timing error is sub-dominant to the achromatic red noise in the post fit residuals. This is because the largest 
contribution to DM variations is a linear trend associated with a DM gradient, which is absorbed into the fitting of 
spin frequency in single-frequency observations. 

A. 2. Determining Observatory and Receiver Backend Offsets 

It is essential to account for phase offsets between telescopes and backend instrumentation that are commonly referred 
to as jumps. Phase offsets were measured relative to the high quality Westerbork ob servations using only overlapping 
data between different data sets and the technique outlined in Uanssen et al.l (|2008l ). Errors in these offsets (possibly 
exacerbated by DM uncertainties) will introduce red noise into the residuals. However, this noise is sub-dominant 
to the observed red noise because it contains both a smaller rms amplitude (as estimated from the uncertainties in 
the jumps) and a shallower fluctuation spectrum because a sequence of imperfectly removed jumps will appear as a 
random walk in pulse phase, which has (approximately) a spectrum cx f~ 2 ■ 

A. 3. Calculating Pulsar Timing Solution 

Using the measured receiver and backend offsets, the combined data set was used to determine the pulsar position 
and proper motion. We did not attempt to measure pulsar parallax. While fitting for astrometric terms, the red noise 
was modeled using higher order frequency derivatives (up to and including d 8 v/dt 8 ). With the astrometric terms fixed, 
a solution including only v and dv jdt was produced that was used in the analysis presented in the main text. 

B. CALCULATING MAXIMUM ENTROPY POWER SPECTRA 

We use a maximum entropy spectral estimator (jPapoulis 198lT) to form power spectra. This type of spectral estimator 
was used because unlike a Fourier transform, it is not affected by spectral leakage. This feature enables us to properly 
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Table 3 

Observations used in this analysis 



Observatory 


MJD Range 


T 


Backend 


V 


^TOA 


°TOA 


Ref. 








(y) 




(GHz) 




(ns) 




Arecibo 


46024 - 


48974 


8.1 


See (1) 


1410, 2380 


420 


200 


1 


Nangay 


47517 - 


53339 


14.9 


VCO 


1410 


3200 


300 


2 


Nangay 


53272 - 


54730 


4.0 


BON 


1370 


194 


30 


3 


EfTelsberg 


51557 - 


54216 


7.3 


EBPP 


1410 


161 


40 


4 


Westerbork 


51566 - 


55375 


10.4 


PUMA 


840, 1380, 2280 


226 


100 


5 


Total 


46024 - 


55375 


25.8 






4201 







Note. — We list the observatory, observation span, backend instrumentation, typi- 
cal observing frequency, typical TOA uncertainty ctoa , and r eferences that describ e the 
instrumental setup. References: (1) Kaspi ct al. (1994) . (2) Cognard ct al. (1993), (3) 
ICognard et~aTI | |2009|) , (4) ILange et al.1 1 12000) . and (5) IVoute et al.l 1 120021)" 



model the red signal observed in the residual TOAs of PSR B1937+21 and produced by the simulations of asteroid 
belts. 

To form the maximum entropy power spectra, we first used cubic spline interpolation to evenly sample the obser- 
vations onto a weekly grid. The maximum entropy spectrum was then constructed using an order 150 autoregres- 
sive model of the interpolated data. A spectrum was calculated using 1000 frequency bins spaced evenly between 
0.038 yr -1 = 1/T and 10 yr -1 . This number of bins provided oversam pling sufficient to detect features in the power 
spectra associated with single dominant periodicities (Press et~all[l99l . 

C. CALCULATING THE TRANSMISSION FUNCTION F(P rb,T) 

For objects with orbital periods P or b much greater than the observing span T, a significant amount, but not all, 
of the signature induced in the TOAs will be absorbed by the fit for pulsar spin down and astrometric terms. We 
quantify the amount of attenuation by a factor F(P OT i,,T). In this section we describe the simulations that we used 
to determine the form of this factor appropriate for our observations. 

We conducted a series of simulated TOAs, with observing cadence matching the actual observations. In each 
simulation, the residual TOAs contained a single sinusoid with period P rb- To account for the TOA modeling process, 
we fit and removed a quadratic polynomial, annual, and semiannual sinusoids from the TOAs. We then calculated the 
ratio of prefit and postfit rms. To produce an estimate of the ensemble average value of _F(-P or b, for each value of 
P or b, we injected the sinusoid a 100 times with different, randomly chosen phases. We calculated this for 200 values of 
P or b for values between 0.2 yr and 210 yr to produce the curve displayed in Figure [3j We used a non-linear fitting 
routine to empirically find the functional form for the transmission function displayed in Equation ([2]). We numerically 
interpolated to produce the single-object sensitivity curves displayed in Figure O 
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